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Observations with intcrfcromctric gravitational- wave detectors result in probability sky maps 
that are multimodal and spread over >10-100 deg^. We present a scheme for maximizing the 
probability of imaging optical counterparts to gravitational- wave transients given limited observ- 
ing resources. Our framework is capable of coordinating many telescopes with different fields 
of view and limiting magnitudes. We present a case study comparing three different planning 
algorithms. We find that, with the network of telescopes that was used in the most recent joint 
LIGO-Virgo science run, a relatively straightforward coordinated approach doubles the detection 
efficiency relative to each telescope observing independently. 

Subject headings: gravitational waves — telescopes — methods: numerical — gamma ray burst: general 



ABSTRACT 



Among the most promising sources for co- 
incident electromagnetic (EM) and gravitation- 
al-wave (GW) emission are compact binary coa- 
lescence (CBC) events in which two neutron stars 
(NSs) or a NS and a black hole (BH) inspiral and 
merge. The NS may be torn apart before merger, 
powering EM emission at a variety of timcscalcs 
and wavelengths ranging from seconds to months 
in 7- and X-ray to radio, respectively (Mcszaros 
2006; Nakar 2007; Mctzger & Berger 2012). It 
is suspected that CBCs are also the progeni- 
tors (Bclczynski ct al. 2006; Leo & Ramircz-Ruiz 
2007a; Troja ct al. 2008) of some or all short, 
hard 7-ray bursts (SHGRBs). Plausible CBC 
event rates suggest that Advanced LIGO (aLIGO; 
Harry & the LIGO Scientific Collaboration 2010) 
and Advanced Virgo (Ad Virgo; Virgo Collaboration 
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Introduction 



2009) could detect about 40 NS-NS and 10 
NS-BH events per year of observation time 
(Abadic ct al. 2010). Aside from shedding light 
on the true progenitors of SHGRBs (Sylvcstre 
2003; Lcc & Ramircz-Ruiz 20071)), coincident de- 
tections of EM and GW emission could also be 
utilized as "standard sirens" to provide a preci- 
sion measurement of the Hubble constant (Schutz 
2001; Nissankc ct al. 2010). 



In addition to CBCs, there are other possi- 
ble sources of coincident EM and GW emission. 
Examples include type II supernovae, which are 
much more common than CBCs at < 1 event per 
20 years in the Milky Way alone, but whose GW 
counterparts would be detectable by Advanced 
LIGO (aLIGO) only within the local group (Ott 
2009). Soft 7 repeaters (loka 2001; Abbott et al. 
2008) and anomalous X-ray pulsars (Abadic ct al. 
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2011) could produce GWs, and some models pre- 
dict GW emission comparable to the best upper 
limits (Corsi & Owen 2011). Neutron stars recov- 
ering from pulsar glitches might produce GW tran- 
sients on timescales of days to weeks that could be 
detectable at <10 kpc (Bennett ct al. 2010). More 
speculative sources of GW include cosmic string 
cusps (Daniour & Vilcnkin 2000; Siemens ct al. 
2006), which if they exist might also power 
GRBs (Babul et al. 1987; Berezinsky et al. 2001; 
Wang et al. 2011) and radio bursts (Vachaspati 
2008). 

The advent of deep, wide-field optical survey 
telescopes and rapid response robotic telescopes is 
paving the way for time-domain astronomy. The 
simultaneous construction of advanced GW detec- 
tors creates the enticing possibility of a coincident 
detection of a GW transient and an associated 
EM counterpart. Efforts in this direction were 
recently made in LIGO's sixth and Virgo's second 
and third science runs, which saw the introduction 
of a program (Kanncr ct al. 2008; Abadic et al. 
2012b, a; Aasi et al. 2012b, a) whereby potential 
GW transient alerts were sent to X-ray, optical, 
and radio facilities for follow-up. In the context 
of space-based (Kocsis ct al. 2008) and advanced 
ground-based (Fairhurst 2009; Cannon ct al. 2012; 
Manzotti & Dietz 2012) GW detectors, there is 
even discussion of using GWs to trigger ear- 
ly-warning follow-up observations while the in- 
spiral is still ongoing. 

The great disparity in the character of optical 
and GW observations makes GW-triggered detec- 
tion of optical counterparts challenging. Although 
there arc now several operating ultra-wide-field, 
but relatively shallow, transient monitoring tele- 
scopes such as Pi of the Sky, only a handful of 
operating or planned deep survey telescopes have 
fields of view (FOVs) that are comparable to the 
error regions of ^10-100 deg^ (Fairhurst 2009; 
Wen t Chen 2010; Nissanke et al. 2011; Fairhurst 
2011) that will be typical of GW sky maps in the 
advanced detector era. By contrast, a GW detec- 
tor is both blessed and cursed by its nature as an 
all-sky instrument: it has coarse directional sensi- 
tivity only due to its geometry-defined quadrupole 
antenna pattern. 

In this paper we lay a framework for plan- 
ning optimal EM follow-up of GW candidates. 
Our paradigm focuses on tiling the GW sky map 



as efficiently as possible by coordinating a num- 
ber of telescopes with dissimilar FOVs. Although 
our formalism supports it, we do not yet account 
for seeing, different limiting magnitudes, different 
times at which observations are made, or multiple 
observations with each telescope. As preliminar- 
ies, we describe the key features of ground-based 
interferometric GW detectors as contrasted with 
optical telescopes. We demonstrate that the opti- 
mal pointing of a single telescope can be expressed 
as a convolution-like integral on the unit sphere, 
which is amenable to spherical harmonic analy- 
sis. For multiple telescopes, we provide a selec- 
tion of observation planning algorithms that take 
into account the particular geometry of each tele- 
scope's image plane as well as local observing con- 
ditions (in our present implementation, just the 
position of the sun and the Earth). We have de- 
veloped a parallelized observation planning code in 
Python/C/C-f— I- called BAYEsian optimal Search 
for Transients with Autonomous and Robotic tele- 
scopes (BAYESTAR)^^, which can map out fol- 
low-up telescope pointings with low latency in the 
way described in this work. We show that with the 
telescopes that were available for EM follow-up in 
the last LIGO- Virgo science run, if each telescope 
performs just one pointing, the imaging efficiency 
is doubled by pointing all of the telescopes in 
coordination rather than independently. Finally, 
we propose an expanded simulation campaign and 
discuss some practical issues for EM follow-up in 
the advanced GW detector era. 

2. GW detectors and optical telescopes 

In this section, we review the observational ca- 
pabilities of GW detectors as contrasted with op- 
tical telescopes. 

GW detectors such as LIGO (Abbott et al. 
2009) and Virgo (Accrncsc ct al. 2008) are Fab- 
ry-Perot Michelson laser interferometers that 



pun on the Cylon battleships in the American television 
series Battlestar Galactica. The defining characteristic of 
the Cylons is that they repeatedly defeat humanity by using 
their superhuman information-gathering ability to coordi- 
nate overwhelming forces. The name also suggests that, 
like the Cylons, robotic transient telescopes may some day 
rise against us humans. 
^We do not like to mention the final 'T' in the acronym, be- 
cause then it would be called BAYES-TART, which would 
sound ridiculous. 



2 



sense GW-induced strain through differences in 
Hght travel time in two arms. The only di- 
reetional sensitivity of a single GW detector is 
due to its antenna pattern (see, for example, 
Sathyaprakash & Schutz 2009). As predicted in 
general relativity, GWs are transverse and come 
in two orthogonal polarizations, the so-called '+' 
and 'x' modes. An interferometer with perpen- 
dicular arms is most sensitive to a polarized 
source at the zenith (or nadir), with the sensitiv- 
ity going to zero at any azimuth that is 45° from 
either arm. The sensitivity to the 'x' mode is 
also greatest at the zenith, but vanishes at any az- 
imuth that is at a right angle to either arm. It also 
vanishes everywhere in the plane of the detector. 

A signal's time delay on arrival at two GW de- 
tectors constrains the position of the GW source 
to an annulus on the sky. The antenna patterns 
of the two detectors impose constraints that dis- 
favor parts of the annulus (Raymoiicl et al. 2009). 
Time delay measurements from three gravita- 
tional wave detectors constrain a source to two 
regions that are located at mirroring points on 
either side of the plane formed by the detectors. 
The mirror degeneracy is, again, partly broken 
by the antenna patterns. Networks of four or 
more detectors are sufficient to determine the 
location of the source to within a single patch 
on the sky. The network consisting of Advanced 
Virgo (AdVirgo) and the two aLIGO observa- 
tories will be able to constrain sources to ar- 
eas of 10-100 deg^ with 90% confidence, varying 
across the sky as a result of the antenna pat- 
terns (Wen & Chen 2010; Nissanke et al. 2011; 
Fairhurst 2011; Khmcnko et al. 2011). Adding 
the planned Japanese detector KAGRA (for- 
merly LCGT; Kuroda 2010) or the proposed 
LIGO-India detector would narrow the confidence 
region (Fairhurst 2011; Schutz 2011). 

Optical telescopes are markedly different from 
gravitational wave detectors in that they have di- 
rectional precision of arcseconds or better. This 
is at the expense of observing only a small frac- 
tion of the sky at one time. Table 1 shows a list 
of telescopes that participated in the LIGO-Virgo 
EM follow-up program. FOVs vary from scarcely 
0.01 deg^ with relatively deep instruments such as 
Liverpool Telescope at a limiting 21 mag; to a few 
deg^ with small robotic transient follow-up tele- 
scopes such as ROTSE at 17 mag; to many deg^ 



with survey instruments such as Skymapper and 
Palomar Transient Factory (PTF), at >20 mag. 
The largest FOVs of >100 deg^ belong to tran- 
sient monitoring instruments such as Pi of the Sky 
at >12 mag. The eagerly awaited Large Synoptic 
Survey Telescope (LSST) will have an even larger 
FOV than PTF, but with a limiting magnitude of 
24.5 in a similar exposure time. Slew times can 
be as brief as seconds for small robotic telescopes 
like ROTSE and TAROT, or minutes and longer 
for larger instruments or instruments that are less 
optimized for response time. Typical total expo- 
sure times range from ~^10 s to a few minutes per 
pointing. 

3. Single telescope case 

In this section, we introduce the optimization 
problem of pointing a single telescope to maximize 
the probability of imaging the true, but unknown, 
location of the GW source. Then, we describe two 
numerical procedures to find the optimum. 

For GW detectors, the observables are strain 
time series consisting of noise plus the astrophysi- 
cal GW signal. The GW signal from a CBC source 
is determined by extrinsic parameters consisting of 
right ascension </>, declination {tt/2 — 9), luminosity 
distance D^, inclination angle t, coalescence time 

polarization angle tp, and orbital phase at coa- 
lescence ifi] and intrinsic parameters including the 
component masses and spins. One may form the 
posterior distribution of all of these parameters: 

p{d, 4>, Dl, t, ip, ^p, mi, TO2, ■ • ■ |gw). 

The proposition GW represents the entire GW ob- 
servation: the strain in each GW detector, the ge- 
ographic locations of the detectors, and their an- 
tenna patterns. If one fixes or marginalizes over 
all parameters except for the position on the sky, 
the spherical polar coordinates u! = (^', ^), then 
this reduces to what we term the GW sky map, 

p{uj\gw) = s{u}). 

This quantity can be evaluated from time de- 
lays and antenna factors (Al)adic et al. 2012a) or 
from the GW strain itself using analytic meth- 
ods (Searle et al. 2008, 2009) or Markov chain 
Monte Carlo (MCMC) techniques (Christensen k Meyer 
1998; van der Sluys et al. 2008; Veitch & Vecchio 
2010; Nissanke et al. 2011). 
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Table 1 

Telescopes Employed in Joint LIGO- Virgo Science Run 



site 


field of view 


limiting 
magnitude 


slew 
time (s) 


geographic location 




LiverpooP*^ 


0.077°x0.077° 


21 


in 


300 s 


30 


28°45'44.8" 


N, 


17°52'45.2" 


W 


Zadko^ 


1.4° X 1.4° 


21 


in 


180 s 


20 


31°21'24" 


s, 


155°42'49" 


E 


ROTSE Ill-a^ 


1.85°xl.85° 


17 


in 


5 s 


4 


31°16'24.1" 


s, 


149°3'40.3" 


E 


ROTSE Ill-b^ 


1.85°xl.85° 


17 


in 


5 s 


4 


23°16'18" 


s, 


16°30'00" E 


ROTSE III-c^ 


1.85°xl.85° 


17 


in 


5 s 


4 


36°49'30" 


N, 


30°20'0" E 




ROTSE Ill-d^ 


1.85°xl.85° 


17 


in 


5 s 


4 


30°40'17.7" 


N, 


104°1'20.1" 


W 


TAROT°f 


1.86°xl.86° 


17 


in 


10 s 


1.5 


43°45'8" 


N, 


6°55'26" E 




TAROT-S°f 


1.86°xl.86° 


17 


in 


10 s 


1.5 


29°15'39" 


s, 


70°43'56" W 


SkyniapperS 


2.373°x2.395° 


21.6 


in 


110 s 




31°16'24" 


s, 


149°3'52" E 




3.5°x2.31° 


20.6 


in 


60 s 




33°21'21" 


N, 


116°51'50" 


W 


QUEST' 


3.6°x4.6° 


21 


in 


140 s 




33°21'21" 


N, 


116°51'50" 


w 


Pi of the Sky SoutW^ 


20°x20° 


12 


in 


10 s 


60 


22°57'12" 


s, 


68°10'48" W 


Pi of the Sky NortW"^ 


40°x40° 


12 


in 


10 s 


40 


37°6'14" 


N, 


6°44'3" W 





Note. — FOV, limiting magnitude, slew time, and geographic location of ground-based robotic optical 
telescopes that participated in the LIGO-Virgo EM follow-up program. Except where otherwise indicated, 
we quote the shortest exposure for which a limiting magnitude is available in the literature. 

^LIGO's request for Liverpool Telescope time calls for a 300 s exposure based on a limiting magnitude of 
21 (White 2012). 

References. — ''Gomboc et al. (2004); ^Guidorzi et al. (2006); ^Coward et al. (2010); ^^Akerlof et al. 
(2003); '^Boer et al. (2003); ^'Klotz et al. (2008); ^Keller et al. (2007); "^Law et al. (2009); 'Baltay et al. 
(2007); JMalek et al. (2010); '^Majcher et al. (2011); 
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Consider the ith of N telescopes labeled 1 . . . A^. 
Let the proposition em,; denote the detection of 
an associated optical transient with telescope i. 
Let Ci represent the configuration of telescope i, 
which consists of a pointing 7^ = {9i,4)i), nominal 
limiting magnitude, exposure start and end times, 
as well as uncontrollable factors such as seeing and 
interference from the horizon, sun, moon, and the 
galactic disk. If wc knew the position u) of the 
source, we could evaluate the probability of de- 
tecting an EM counterpart at u) with telescope i, 

p{¥Mi\Ci,Uj). 

For most of this paper we shall restrict ourselves 
to the dependence on the pointing of the telescope, 
so that the detection probability factors into a 
function hi{ijj) describing the telescope's FOV in a 
home orientation (i.e., for equatorial mounts, the 
celestial North pole) and a function Viiuj) that de- 
scribes environmental factors including the night- 
time/twilight terminator and the horizon. In prin- 
cipal, Vi{ijj) could also include other observational 
factors such as moonlight and current seeing. The 
detection probability becomes: 

p(EM,|7,,w) {R~\'^i]u:)v,{u}). (1) 

A simple rectangular FOV with angular width w 
and height h is described by 



if I sin6'cos(/i| < sinw/2, 
|sin0sin(/)| <sin/i/2, 
and 6 < 7r/2; 
otherwise. 



(2) 



Much more complicated FOVs are commonplace 
because a real telescope may have a focal plane 
consisting of a CCD mosaic with a filling factor 
that is much less than unity; dead CCDs or CCDs 
that contain defects; or CCDs that are vignetted 
or clipped by the optics. 

Rl'Ji] denotes an active rotation operation 
about the z-axis by 0^ and about the y-axis by 
6i, so -R~^[7i] is the passive rotation such that 
bi evaluated at (i?r^[7ja;) represents the FOV 
centered on the sky location 7^ = {Oi,(f)i). This 
rotation matrix is described by two, not three, an- 
gles because mounts for ground-based telescopes 
usually have two degrees of freedom (for example, 
right ascension and declination or azimuth and 
altitude) . 



This particular model incorporates the assump- 
tion that a telescope always detects an optical 
transient when the true location of the GW source 
is inside the FOV. We can identify this with the 
probability of imaging the source, which for most 
of this work we shall use as a proxy for the ulti- 
mately desirable probability of detection. 

The optimal configuration C* of telescope i 
maximizes the probability of detection of an elec- 
tromagnetic counterpart given all gravitational 
wave observations, marginalized over the unknown 
source location. To wit, 



C* = are max 



p(EM,|C,,a;)p(u;|Gw)df2, (3) 



where the integral is over uj and AO, represents the 
surface element on the unit sphere. In the simpli- 
fied scenario in which the pointing orientation is 
the only free parameter, the optimal pointing is 
given by 

7* = argmax j p(EMi|7j, a;)p(a;|Gw) dfJ 

= argmax j s(a;) 5; (i?^^[7,j] a;) 11^(0;) dfl. 

(4) 

The integral in Equation (4) resembles a convo- 
lution integral, but over rotations instead of over 
the real number line. Just as fast Fourier trans- 
form (FFT) methods are used to efficiently evalu- 
ate convolutions, fast spherical harmonic methods 
may be used to evaluate Equation (4) efficiently. 
In the Appendix, we describe in detail two algo- 
rithms for solving the single telescope problem: a 
"direct" or spatial algorithm and a "multipole" 
algorithm based on the method for fast convolu- 
tion on the sphere of Wandclt & Gorski (2001). 
Both algorithms exploit special properties of the 
HEALPix^ (Gorski et al. 2005) projection. Either 
of these algorithms may be used as a component 
of the planning of observations with multiple tele- 
scopes, which we discuss in the next section. 

4. Multiple telescope case 

In order to address the problem of pointing mul- 
tiple telescopes, we must first specify a figure of 



^http : //healpix . jpl . nasa . gov/ 
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merit. One example is the probability of observ- 
ing an electromagnetic counterpart in any (one or 
more) telescope, 

Pany = ^(eMiUEMsU- • -UEMjv |7l , 72 . • ■ - 7Ar, Gw). 

(5) 

This only considers the case in which each tele- 
scope is asked to point once. Some telescopes will 
take a sequence of images at arbitrary sky loca- 
tions. Others are constrained to follow a particu- 
lar path on the sky, for instance, visiting reference 
pointings in a specified order. In the former case, 
we could trivially interpret eMi . . . EM^v as opti- 
cal detections in the observations 1 ... at times 
ti . . An respectively, where there are N observa- 
tions but potentially ^ TV telescopes. In the lat- 
ter case, we could define a composite FOV that 
describes the probability of imaging the source in 
any of the pointings throughout the entire obser- 
vation. A more elaborate objective function would 
also include a penalty on the number of pointings 
in order to economize the limited observational re- 
sources. For the remainder of this paper, we will 

focus on Pany**- 

De Morgan's law tells us that Equation (5) is 
equivalent to the complement of the probability of 



* Alternatively, we could demand not one, but several, EM 
detections to secure the discovery of a bona fide optical 
transient. Combinatorial considerations give the probabil- 
ity of detection in at least K out of A'^ observations, 

j)(at least K detections | Gw) = 

L i 
+ („l)^f+2 (^^^ ^ J Ep(EMi n EMj |7i, 7i, 



A - In 
K -l) 



P{eMi n EMiv|71, 



p(a>|Gw) df2. 



,7iv,tj) 
(6) 



non-detection with every telescope, 

Pany = 1- j [1 ^ P(EM,|7i, w)]^ p(u;|GW) dfj 

• • • [1 - 6jv {R'^ilN] ^) vnM] s{uj) dn. 

(7) 

The optimal set of pointings maximizes Pany, 

(7l,72>--->7w) = argmax pany (8) 

7i,72.---7jv 

4.1. Planning algorithms 

We considered three different multiple telescope 
algorithms, increasing in sophistication and accu- 
racy but also increasing in computational over- 
head. After describing each algorithm in detail, 
we exhibit implementations that use the positional 
astronomy library NOVAS^ (Barron et al. 2011) 
to take into account rudimentary sun and horizon 
interference assuming observers at sea level on an 
atmosphereless, perfectly spherical Earth. We as- 
sume that all observations occur simultaneously at 
the time of the GW detection. 

4-. 1.1. Noncoordinated 

The noncoordinated planner models the situ- 
ation in which each telescope operator indepen- 
dently decides to observe that patch of the sky 
that maximizes his or her own odds of detecting 
a counterpart. Each single telescope problem is 
solved independently, regardless of the configura- 
tions of any of the other telescopes. 

The noncoordinated method has the disadvan- 
tage that it will generally result in all of the tele- 
scopes pointing at the same region of the sky, per- 
haps failing to image all of the significant regions 
of the GW sky map. However, if the telescope net- 
work's coverage is exceedingly poor, such that the 
nonzero regions of the sky map are within view of 
only ~ 1 telescope at a time, then this will not be a 
significant drawback. The noncoordinated method 



http; / /aa.usno .navy .mil /software/no vas/novas_ info .php 
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also has the advantage that it is exceedingly sim- 
ple to implement. For these reasons, the noncoor- 
dinated planner may represent a reasonable, but 
sub-optimal, observational strategy. 

4.1.2. Greedy 

The greedy planner is slightly more sophisti- 
cated than the noncoordinated method. Suppose 
that pointings have been decided upon for tele- 
scopes 1,2, ... ,i. Then the configuration of tele- 
scope i 1 is chosen such that it maximizes the 
imaging probability when all of the previous point- 
ings 1 . . . i are held fixed. This procedure is illus- 
trated in Figure f. It is not guaranteed to find 
an absolute maximum, but we shall see that in 
practice it does quite well. 

This algorithm leaves a choice of the order in 
which to point the telescopes. This is analogous 
to the problem of packing a number of stones of 
different sizes into a bucket. The intuitive an- 
swer is to insert the largest rocks first and then 
use the pebbles to fill in the gaps. Likewise, we 
find that we obtain the best detection efficiency 
when we sort the telescopes by decreasing FOV 
area because the telescope with the largest FOV 
makes the greatest progress toward imaging the 
entire support of the GW sky map. We refer 
to the greedy algorithm with this prescription as 
the greedy-sorted planner. We expect that the 
greedy-sorted planner will always outperform the 
noncoordinated planner by a factor that depends 
on the signal-to-noise ratio (SNR) of the GW sig- 
nal and the particular network of telescopes. 

4-. 1.3. Simulated annealing 

The anneal planner uses the simulated anneal- 
ing code provided by the SciPy library^ to ex- 
tremize Pany Starting from an initially proposed 
observing plan, this procedure randomly perturbs 
all of the telescopes' pointings simultaneously. If 
the imaging probability has improved, the per- 
turbed plan is adopted. The perturbations be- 
come smaller and smaller with each iteration; the 
perturbation distribution as a function of time is 
called a cooling schedule. We used a modified 
version of the "very fast" cooling schedule de- 
scribed by Ingber (1989) with an initial "temper- 



'http : //docs . scipy . org/doc/scipy/ref erence/generated/scipy . 



( \ 

Begin 



Read next 
telescope's FOV 



Mask parts of sky map 
that are in sunlight 
or blocked by Earth 




yes 



Convolve FOV 
with sky map 




Mask out FOV 
from sky map 




no 



[Pone] 

ipfigizi-.ann&imiffiliart depicting greedy algorithm. 
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ature" Tq = 10~^. Since Ingbcr's cooling schedule 
uses rectilinear coordinates, we modified it by per- 
turbing the Cartesian coordinates of the center of 
each telescope's FOV and then normalizing to unit 
length. 

Simulated annealing will nearly always find the 
globally optimal observing plan in a finite, but 
possibly large, number of iterations. With a par- 
ticular fixed number of iterations, we expect that 
it will usually outperform not only the noncoordi- 
nated planner but also the greedy-sorted planner. 
However, the computational overhead of evaluat- 
ing thousands of observing plans in sequence may 
be prohibitive if there is a time constraint — if, for 
example, it is important to begin taking images as 
soon as possible after the GW event. 

4.2. Case study 

We compared the imaging efficiency of our three 
multiple telescope planners using a set of 2126 GW 
sky maps from signals injected into simulated ini- 
tial LIGO and Virgo noise over a period of 24 
hours. The simulated GW signals are accurate 
to second post-Newtonian order in phase. Compo- 
nent masses arc drawn from a uniform distribution 
from 1 to 15 Mq with total mass is restricted to 
<20 AIq. Distances are drawn randomly from a 
log-uniform distribution from 10 to 40 Mpc. A 
sky map was produced for each injection using 
the time delay method described in Abadic et al. 
(2012a). We generated two observing plans for 
each sky map for each of the three planners. The 
first observing plan used Vi{u:) = 1 for all observa- 
tions i and sky locations u;, disregarding the posi- 
tions of the sun and Earth. The second observing 
plan used Vi(u}) = 1 above the horizon and at 
all sky locations from which the sun was at least 
18° below the horizon, and zero elsewhere, simu- 
lating observing conditions on an atmospherelcss, 
spherical Earth. The output of the greedy-sorted 
planner was used as the initial proposal for simu- 
lated annealing. Pi of the Sky and Liverpool Tele- 
scope were excluded from the simulation because 
the areas of their FOVs are outliers with regard 
to the other FOVs of ~1~10 deg^ represented in 
the model. Figure 4 illustrates an observing plan 
generated with the greedy algorithm with sun and 
horizon constraints disabled. 

Figure 2 shows the measured imaging efficiency 
as a function of the luminosity distance. For each 



sky map and for each planner, we tested whether 
the injected location would have been imaged, i.e., 
whether the injected source location was inside 
any of the FOVs of the telescopes. Dividing the 
injection set into six concentric shells, each 5 Mpc 
deep, from 10 to 40 Mpc, we tallied the number 
of sources that would have been imaged and the 
total number of sources. Because the injections 
were drawn from a uniform distribution in log D l , 
we weighted each injection by 

-1 



Dr 



d log Dl 



(9) 



in order to approximate a uniform-in-volume dis- 
tribution. The dashed lines represent the observ- 
ing plans that disregard the presence of the sun 
and the Earth whereas the solid lines represent 
the observing plans that fully account for these 
effects. 

For each injection, we also computed the proba- 
bility Pany of imaging the source, marginalized over 
(6*, (/)) and not making use of knowledge of the true 
sky location, using Equation (7). We then formed 
a histogram of pany over all of our injections, again 
weighting by the correction factor D^. This his- 
togram is shown in Figure 3. 

Two systematic effects are apparent in Figure 2. 
First, for all planning algorithms, accounting for 
the sun results in a reduction in efficiency by a fac- 
tor of sa 1/3, independent of distance. This is be- 
cause astronomical twilight begins and ends when 
the sun is 18° below the horizon (Astronomical Al- 
manac 2010); solar illumination is faint enough for 
astronomical observations from the Earth's sur- 
face over over a fraction (1 — cos [90° — 18°]) /2 w 
0.35 of the sky. Second, the efficiency falls as lu- 
minosity distance increases because the area of the 
GW error region grows as signal-to-noise ratio de- 
creases. 

Both Figures 2 and 3 demonstrate that co- 
ordination of many optical telescopes increases 
the chances of imaging the true location of the 
GW source. Whether or not sun and hori- 
zon interference are included, the anneal and 
greedy-sorted planners achieve about twice the 
detection efficiency of the noncoordinated planner 
for the network of telescopes used in the last joint 
LIGO-Virgo science run and in our model. 

From Figure 3, we see that simulated anneal- 
ing generally outperforms the greedy algorithm. 
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Fig. 2. — Average detection efficiency of nonco- 
ordinated, greedy-sorted, and anneal planners, for 
sources located within spherical shells 5 Mpc thick 
at 10-40 Mpc. Solid lines include sun and horizon 
interference; dashed lines ignore interference con- 
siderations. 
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Fig. 3. — Histogram of expected detection effi- 
ciency Pany as coniputcd using Equation (7), for 
sources uniformly distributed in volume from 10 
to 40 Mpc. As in Figure 2, solid lines include sun 
and horizon interference and dotted lines ignore 
interference considerations. 



At least a slight improvement is expected because 
the observing plan generated by the greedy-sorted 
algorithm is used as the initial proposal for simu- 
lated annealing. Simulated annealing would also 
give only a slight improvement if the greedy-sorted 
observing plan was already nearly optimal. 

5. Discussion 

We have presented a framework for planning 
optimal EM follow-up of GW transients. We have 
shown that, if each telescope images one sky lo- 
cation, a coordinated approach with a particular 
network of telescopes is twice as likely to image the 
true location of a GW event, versus an uncoordi- 
nated approach. It remains to be seen whether 
coordination will be important for more realis- 
tic observing strategies, which will involve tak- 
ing multiple images with each telescope. A suc- 
cessful observational program also will likely in- 
volve several different stages calling for markedly 
different kinds of instruments. In the context of 
a grand strategy for the follow-up of GW candi- 
dates (Metzger & Berger 2012), spatially coordi- 
nated observing campaigns like those we have con- 
sidered are most likely to be applicable to prompt, 
rapidly fading, on-axis afterglows, though it is 
unclear whether coordination will be relevant for 
later-rising and fainter counterparts. This will be 
the subject of future work. 

We have found that a simple greedy algorithm 
performs almost as well as simulated annealing 
while taking much less time on average. In the 
interest of low latency, we have provided a par- 
allel implementation of the convolution algorithm 
so that observation plans can be generated rapidly 
on a multicore workstation. The observation plan- 
ning code has been operated in a low latency 
mode during software engineering runs of the joint 
LIGO-Virgo low latency data analysis. 

In this work, the objective function that we 
have used is the probability of imaging the source. 
This is a proxy for the more desirable objective 
function of the probability of detecting the source, 
which must take into account the limiting mag- 
nitudes of the observations, a description of the 
potential false positives, and a model of the EM 
emission. A weighting could be assigned to each 
telescope's FOV based upon the integrated flux 
that the telescope would be able to gather given 
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the anticipated power law index of the light curve 
and the telescope's limiting magnitude, slew time, 
and exposure time. For example, assuming a pow- 
er-law light curve i^^ oc i^", at time t a telescope 
with limiting magnitude M would be able to de- 
tect a source up to a maxim um luminosity dis- 
tance Dl max OC V2MWr^. We can marginalize 
over the unknown distance to the source, assum- 
ing that nearby 7-ray bursts (GRBs) are uniformly 
distributed in space. Neglecting cosmological cor- 
rections, the detection probability of a single tele- 
scope might be expressed by 



p(eMj|7,,u;) cx 




")^^^ inside FOV, 
outside FOV, 
(10) 

where ti is the time of observation i and Mi is 
the limiting magnitude of observation i. We could 
also marginalize over the power-law index a and 
the initial luminosity. 

Although this prescription will be appropriate 
for a single telescope, for multiple telescopes or 
multiple pointings of the same telescope it will be 
more appropriate to form the probability of detec- 
tion with the network of telescopes, conditioned 
on the luminosity distance, and then marginalize 
over luminosity distance. Cast in this form, the 
multiple telescope detection probability is 



Some of the results in this paper have been de- 
rived using HEALPix (Gorski et al. 2005). 

LIGO was constructed by the California Insti- 
tute of Technology and Massachusetts Institute of 
Technology with funding from the National Sci- 
ence Foundation (NSF) and operates under co- 
operative agreement PHY-0107417. Some results 
were produced on the NEMO computing cluster 
operated by the Center for Gravitation and Cos- 
mology at University of Wisconsin-Milwaukee un- 
der NSF Grants PHY-0923409 and PHY-0600953. 
This research is supported by the NSF through 
a Graduate Research Fellowship to L.S. A.S. ac- 
knowledges support from the LIGO Laboratory 
Summer Undergraduate Research Fellowship, Vis- 
itors, and NSF Research Experiences for Under- 
graduates (REU) programs. This paper has LIGO 
Document Number LIGO-P 1200033- v4. 



Pany 



= 1-- 



1-e [DL,max,» {L^, M,)-Dl] h [R-^h,] w) V^{u:) 



s(a;) 



p(Li, L2, ■■■,Ln) dLidLz . . . dL^ Dl^ADl dO 

(11) 

where V is the volume of integration and O is the 
Heaviside step function. D^^^s^^i cx ^^2.512*^'^^ 
is the farthest detectable luminosity distance for 
observation i, being a function of the source's lu- 
minosity Li at the time t; of the observation and 
the limiting magnitude Mi. This figure of merit 
is not as simple as Equation (7), but it may be- 
have better under simulated annealing because the 
FOVs no longer have to move rigidly around each 
other without overlapping. This is also the subject 
of future work. 



Source code for BAYESTAR is available on the 
LIGO Data Analysis Software Working Group web 

site at http : //www. Isc-group .phys .uwm. edu/daswg/projects/bayestar . html. 
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Fig. 4. — Demonstration of a BAYESTAR observation plan for a simulated NS-NS CBC event injected into 
synthetic initial LIGO noise. The event produced triggers in the HI, LI, and VI detectors. The GW sky 
map, p{u)\gw), in posterior probability per square degree, is shown in color in the online version and in 
shades of gray in the print version. The FOVs of the telescopes are shown as shaded gray rectangles. Panels 
(a) and (b) show Lambert azimuthal equal area projections of two 30° x 30° patches where observations were 
planned. Panel (c) provides the selected telescope pointings in tabular form. Panel (d) shows the selected 
observations on an all-sky MoUweide projection. 
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A. Implementation 

Equation (4) is reminiscent of a convolution integral. This leads us to think of bi{u}) as an analogue of a 
finite impulse response (FIR) filter, but on the unit sphere. A FIR filter with a sufficiently compact kernel 
can be efficiently evaluated in so-called direct form by advancing the kernel one sample point at a time and 
taking an inner product with the signal. On the other hand, a FIR filter with a long, extended kernel is 
generally more efficiently implemented using frequency domain techniques such as the classic overlap-save 
algorithm (Numerical Recipes, Press et al. 2007), which exploits the convolution theorem and the FFT. 

Similarly, if the telescope's FOV spans sufficiently little area, then a discrctized form of Equation (4) can 
be efficiently evaluated by subdividing both the FOV and the GW sky map into pixels, and for each pointing 
•Yj, rotating the pixels in the kernel and then taking their inner product with the pixels in the GW sky 
map. Transforming just the nonzero pixels of the kernel results in a substantial computational speedup; for 
a ~1 deg^ FOV, this reduces the number of calculations by a factor 4(180)^/7r ~ 4 x 10^. 

If, however, the telescope's FOV subtends a large solid angle, then the "direct" algorithm outlined above 
will be computationally expensive, with a cost growing as ©(npix^), where ripix is the number of pixels cover- 
ing the entire unit sphere. The FFT-based frequency domain techniques for FIR filters also have analogous 
fast multipole transform algorithms for functions defined on S^. Considerable study has been devoted to 
fast multipole methods because of their utility in solving differential equations with approximate spherical 
symmetry. Example applications can be found in fields as varied as cosmic microwave background (CMB) 
mapping (Gorski et al. 2005), numerical astrophysics, quantum chemistry, and even the entertainment in- 
dustry (McEwen & Wiaux 2011). Presently, the best known upper bound on the complexity of both spher- 
ical harmonic transforms and fast convolution on the sphere (Gorski et al. 2005; Wandelt & Gorski 2001) 
is ©(ripix^^^) ^ 0{L^), where L ~ 0(npix^^^) is the azimuthal index of the highest non-vanishing (or 
non- negligible) spherical harmonic component. 

For a sufficiently dense pixel resolution and sufficiently large FOV, the multipole method will inevitably 
be more efficient than the direct method. Anticipating that one approach or the other might be more 
computationally efficient for a given FOV, we implemented both a "direct" algorithm and a "multipole" 
algorithm, described below. 

A.l. "Direct" algorithm 

For the "direct" algorithm, it is important that the sky is divided into patches, or pixels, of equal area. 
The pixelization scheme defines a surjective map from points on the unit sphere to integer pixel indices, 
/: !->■ Z„pj^. It is desirable to be able to compute the discretized version of the rotated FOV, bi{'j~^u}), 
by transforming the pixel coordinates themselves. If some pixels are much larger or smaller than others, then 
some rotations will require a great deal of sub-pixel interpolation. An example of a particularly troublesome 
pixelization scheme would be a division into a equiangular lattice in right ascension and declination, because 
pixels near the poles would subtend much smaller solid angles than pixels near the equator. It is much better 
if the pixels are of equal or nearly equal size, for then a "nearest neighbor" interpolation scheme consists of 
nothing more than a permutation of the pixel indices. Such an interpolation scheme ought to converge to 
first order in the pixel area, while being exceptionally fast. 

For this reason, we chose to build an implementation of the "direct" algorithm in HEALPix (Gorski et al. 
2005), a special pixelization of the sphere and an accompanying software package that has become one of the 
workhorses of all-sky astronomical data analysis. HEALPix stands for Hierarchical Equal Area isoLatitude 
Pixelization. As the name says, it has three chief virtues. The first is that successive levels of resolution 
are derived by subdividing pixels. This makes it straightforward to downsample (decreasing resolution) or 
upsample (increase resolution) a map. Also, for search operations it provides a log-linear speedup that is 
analogous to binary searches trees or octrees. The second virtue is that all HEALPix pixels at a given 
resolution have the same area, satisfying the requirement for our "direct" convolution algorithm. Finally, 
the isolatitude property, that pixels are arranged on rings of constant latitude, facilitates computationally 
efficient multipole transforms at any resolution, which we exploit in the next section. 
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Algorithm 1 "Direct" convolution for single telescope case. 

{si[0], Si[l], . . . ,Si[npix — 1]} pixelized map 
{vi[Q\, . . . , Vi[n-p\y^ — 1]} pixelized mask 
{6i[0], ■ • ■ , &i[«pix — 1]} pixelized FOV 
{PiPiMl], . . . ,K["pix - 1]} ^ {0, 0, ... ,0} 
for j = — ripix — 1 do 
7 ^ pix2ang (j) 

for /c G Z: < fc < Tipi^, bi[k] ^ do 

U3 ^ pix2ang (fc) 

k' ang2pix(i?[7-i]u;) 

Pi[j] ^ Pi[j] + bi[k] X Vi[k'] X Si[fc'] 
end for 
end for 



A. 2. "Multipole" algorithm 

The "multipole" algorithm is the analog of FFT convolution in spherical coordinates. Our implementation 
of fast convolution with HEALPix is based on Wandclt & Gorski (2001). A weighted outer product is formed 
from the coefhcients of the spherical harmonic expansions of the sky map and the telescope FOV. The 
resultant tensor is the 2D discrete Fourier transform (DFT) of the convolution on a equiangular grid in 
{0,(j)). Our contribution is to specify a computationally efficient procedure for obtaining the convolution 
on the HEALPix grid instead of the equiangular grid in spherical coordinates. This is an especially useful 
feature for our application because output in HEALPix coordinates facilitates further spatial manipulation 
after the convolution. A detailed explanation of the algorithm is provided below. Source code in C++ is 
provided in the supplementary material. 

Without loss of generality, from here on we drop the environmental factor Vi(uj) and the telescope index 
z, leaving just 

p(em|7,Gw) = J s(u})b [R^^[y]u}) dfl. 

Though the integral involves a passive rotation of the beam by R~^['y], it is conceptually simpler to think of 
this is as an active rotation by Rl'f]. In the language of spherical harmonics, this integral can be evaluated 
at any 7 = {9, (f) through the sum 

L I I 

p(EM|0,<^,GW)=^ ^ S*i„,DlM,0,mn (Al) 

rn— — l n— — l 

where 

DLnia, /3, 7) = e-™-dL„(/3)e-™^ (A2) 

is the Wigner matrix element denoting the active rotation first by 7 about the z-axis, then by (3 about the 
new y-axis, and finally about a about the final z-axis. (Sakurai & Tuan (1994) use the active convention; 
Edmonds (1957) and Mathematica employ the passive convention.) sim and bin arc the spherical harmonic 
coefficients of s{u:) and b{uj), respectively, assumed to be zero or at least negligibly small for / > L. (i^„(/3) is 
the reduced Wigner matrix element representing an active rotation about the y-axis by /3. The basic idea is to 
employ a Fourier expansion of the reduced Wigner coefficients (Edmonds 1957; Risbo 1996; Wandclt & Gorski 
2001), 

I 

4„(/3) = *"-" E AL.™A^,„e-™^ (A3) 

m'——l 
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where Aj„„ = (iJ„„(7r/2), to convert the summation in Equation (Al) to a two-dimensional, discrete Fourier 
transform, 



e ^Oln 



p(EM|0,,/>,Gw) = ^ ^ ^ ^ ^™-"sLAU.AL'„e- 

/— m— — l m' — ~l n— — l 
L L 

^ Tmm'e-"'"'''e-''^* (A4) 

m— — L ni' — ~L 

First, compute both sim and bim using the HEALPix routine map2alm for m > 0. Since 5(0;) and b{u:) 
are purely real, the remaining coefficients to < are determined by conjugate symmetry, si^-m = (^l)™*;*™- 
Also, all of the m = coefficients are purely real. 

Next, compute sim = {—iY^sim and bim = {—i)™'bim for to > to transfer the phase factor i™^" to the 
beam and sky map. 

Next, evaluate the sum over n in 0(L^) time, yielding Xim'- 



Y - \^ Al h - 7, 4-9 A' /R-C M if - m - n) IS even 
^ ^ [Im 6;„ if (/-TO-n) isodd 



(A5) 



The expression after the equality sign makes use of 6;„ ~ and the identity (McEwen & Wiaux 2011) 

dLni^ -fi) = J/3) =^ A:„„ = (-1)'-™A5„^_„.' 

Next, for < TO < L and —L < m' < L, the sum over I is computed, also in 0{L^) time: 

T =V"«* A' Y _^r«* A' /^'"'' ifm'>0 

/— max(|m|,|m'|) /— max(m, |m' | ) v ^l^'^ 1 

At this point, it is possible to recover the convolution on a equiangular grid on right ascension and 
declination by performing a 2D forward Fourier transform of Tmm'- However, this grid will oversample the 
poles; it is preferable to obtain the convolution at the HEALPix pixel sites. The HEALPix grid is arranged 
such that the pixel centers {9j,(j)jk) lie on a set of 0{L) iso(co)latitude rings at colatitudes 9j. For the 
colatitude 9j particular to each ring, the sum over to' is performed in a total of 0{L'^) operations, yielding 

L 

T,n{0j)= J2 T,nm'e-"-^''^. (A7) 

m' — — L 

This ID Fourier series, Tm{Oj), is now shifted in phase so that To(6'j) is the value at the first pixel in the 
ring centered at {(j)jQ^6j), resulting in 

TUe,)^Tmi0,)e-^^*'°. (A8) 

Since the convolution must be everywhere real, we know that T'_^-^{9j) = [27„(0j)]*. 

Now, in rings close to the poles, the number of pixels per ring ?T.ring,j is small; we assume that in all rings 
'T-ring.i < 2L + 1. For more efficient evaluation, the Fourier series can be folded, or aliased, to T," (6*^), for 
< TO < nringj/2, such that 



27rifcml ""B£/ ^.,,„ , \ 2'Kikm 



m=-(nri„gj/2-l) 



''ring,j 



'Ting,:; 



(A9) 
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where 

TM) = E (^r'n+„„„,.,, + 7^™-„.„,„.) for < m < 7iH„g,,/2. (AlO) 

i 

The Fourier series can be summed directly in a total of 0{L^) operations for all rings. Although not 
necessary to preserve the overall cubic scaling with L, this final step can be sped up by employing a FFT 
instead of a direct sum. FFT packages such as fftw^ conventionally provide a real-to- complex forward FFT for 
real-valued input data, producing just the nonnegative frequency components, and a corresponding inverse 
complex-to-real FFT converting just nonnegative frequency components back to real data. We would like 
to use a complex-to-real transform to fill in the ring pixels from the Fourier series, but the negative sign 
in the exponent of Equation (A9) makes it a forward, rather than an inverse, transform. Application of 
a complex-to-real FFT would result in the ring pixels coming out in reversed order. For this reason, we 
compute 

t:::{0,) = ku^,)]* (ah) 

and then perform the complex-to-real FFT on this to fill in the pixels in each ring. 
A. 3. Comparison 

We implemented both the "direct" and "multipolc" convolution algorithms in C/C++. The outermost 
loops are accelerated using OpenMP^, a widely supported compiler extension for generating parallel code for 
multiprocessor machines. Much of the HEALPix library itself, including the spherical harmonic transform 
routine map2alm, is also accelerated with OpenMP. For the multipole algorithm, the final DFTs are performed 
using fftw. 

Operating with one thread on a quad core Sun Fire X2200 M2 and a pixel resolution of «0.05 deg^, the 
direct algorithm ran in 8.9 s for the ROTSE III FOV, 38.4 s for QUEST, and 931.2 s for Pi of the Sky 
South. The multipole algorithm, whose run time does not depend on the FOV, took 23.4 s. Using 8 threads, 
ROTSE-III took 1.3 s, QUEST 5.4 s. Pi of the Sky South 128.9 s, and the multipole algorithm took 6.6 s. 



^http : //www . f f tw . org/ 
*http : / /opermip . org/ 
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